• Maize data taken from a Genome Research study (Paschold et al. 2012).
• 4 genetic lines available:
  • B73 (parent1)
  • Mo17 (parent2)
  • B73xMo17 (hybrid)
  • Mo17xB73 (alternate hybrid, discarded)

\[\begin{align*} \\ &\text{HPH gene $g$ if} \ &&\mu_{g, \ \text{hybrid}} > \max \left ( \mu_{g, \ \text{parent 1}}, \ \mu_{g, \ \text{parent 2}}\right ) \\ &\text{LPH if} &&\mu_{g, \ \text{hybrid}} < \min \left ( \mu_{g, \ \text{parent 1}}, \ \mu_{g, \ \text{parent 2}}\right ) \\ \\ \text{Equ}&\text{ivalently,} \\ \\ &\text{HPH if} &&\alpha_g + \delta_g > \ \ \ |\alpha_g - \delta_g| \\ &\text{LPH if} &&\alpha_g + \delta_g < -|\alpha_g - \delta_g| \\ \end{align*}\]

• Data \(\boldsymbol{y}\), parameters \(\boldsymbol{\theta}\).

• Sample \(\boldsymbol{\theta}\) from \(p(\boldsymbol{\theta} \ | \ \boldsymbol{y})\) by iterating over the Gibbs steps.

  • Sample \(\nu_\rho\) from its full conditional density \(p(\nu_\rho \ | \ \cdots)\) with slice sampling.
  • \(\nu_\gamma\) similar
  • \(\varepsilon_{n, g}\)'s similar
  • \(\rho_n\)'s similar
  • \(\phi_g\)'s similar
  • \(\alpha_g\)'s similar
  • \(\delta_g\)'s similar

• Sample \(\tau_\rho^2\) from its full conditional gamma density with slice sampling.
• \(\gamma_g^2\) similar, but from an inverse-gamma
• \(\sigma_\phi^2\) similar, but from a truncated inverse-gamma
• \(\sigma_\alpha^2\) similar
• \(\sigma_\delta^2\) similar
• Sample \(\theta_\phi\) directly from its full conditional normal density.
• \(\theta_\alpha\) similar
• \(\theta_\delta\) similar

• Run a pilot MCMC chain for a fixed number of iterations.
• Initialize 3 additional chains using dispersed starting values relative to the pilot chain MCMC samples.
• Continue the 4 chains independently until all Gelman-Rubin potential scale reduction factors for "parameters of interest" are less than 1.1 (Gelman, Andrew et al. (2013), Gelman and Shirley (2011)).
  • "Parameters of interest": hyperparameters, \(\alpha_g\)'s, \(\delta_g\)'s.
• After apparent convergence, combine the four chains and resume the MCMC until 100 effective samples are obtained (Gelman, Andrew et al. 2013) for each hyperparameter.

• \(M\) Monte Carlo iterations, index \(m\).
• \(\alpha_g^{(m)}\) and \(\delta_g^{(m)}\) are Monte Carlo samples.
• \[ \begin{align*} p\left( \text{$g$ is HPH} \ | \ \boldsymbol{y} \right) &\approx \frac{1}{M} \sum_{m = 1}^M I\left (\alpha_g^{(m)} + \delta_g^{(m)} > \ \ \ \left |\alpha_g^{(m)} - \delta_g^{(m)} \right | \right) \\ p\left( \text{$g$ is LPH} \ | \ \boldsymbol{y} \right) &\approx \frac{1}{M} \sum_{m = 1}^M I\left (\alpha_g^{(m)} + \delta_g^{(m)} < -\left |\alpha_g^{(m)} - \delta_g^{(m)} \right | \right) \\ \end{align*} \]

• \[ \begin{align*} \widehat{\alpha}_g = \frac{1}{M} \sum_{m = 1}^M \alpha_g^{(m)} \qquad \qquad \widehat{\delta}_g = \frac{1}{M} \sum_{m = 1}^M \delta_g^{(m)} \\ \\ \end{align*} \]
• \[ \begin{align*} \text{Effect size}_g &\approx \left (\widehat{\alpha}_g + \widehat{\delta}_g - \left |\widehat{\alpha}_g - \widehat{\delta}_g \right | \right)^+ \\ &+ \left (\widehat{\alpha}_g + \widehat{\delta}_g + \left |\widehat{\alpha}_g - \widehat{\delta}_g \right | \right)^- \end{align*} \]

1. Embarrassingly parallel Gibbs steps for gene-specific parameters.
2. Parallel reductions for the rest.

• \[ \begin{align*} \\ p(\gamma_g^2 | \cdots ) &= \text{Inverse-Gamma} \left ( \frac{N + \nu_\gamma}{2}, \ \frac{1}{2} \left (\nu_\gamma + \sum_{n =1}^N \frac{\varepsilon_{n, g}^2}{\rho_n^2} \right ) \right ) \\ \\ \end{align*} \]

• \(G = 39656\) parallel GPU threads, one for each \(\gamma_g^2\).

• \[ \begin{align*} p(\varepsilon_{n, g} \ | \ \cdots) &\propto \exp \left ( \varepsilon_{n, g}y_{n, g} - \frac{\varepsilon_{n, g}^2}{2 \rho_n^2 \gamma_g^2} - \exp(\varepsilon_{n, g}) \exp(\eta(n, g)) \right) \\ \\ \end{align*} \]

• \[ \begin{align*} A &= y_{n, g} \\ B &= \frac{1}{2 \rho_n^2 \gamma_g^2} \\ C &= 0 \\ D &= \exp(\eta(n, g)) \\ E &= 0 \\ \\ \\ \end{align*} \]

• \(N \times G = 12 \times 39656\) threads, one for each \(\varepsilon_{n, g}\).

• \[ \begin{align*} &\frac{d^2}{d \nu_\gamma^2} \log p(\nu_\gamma \mid \cdots) = - \frac{G}{4} \psi^{(1)} \left ( \frac{\nu_\gamma}{2} \right) + \frac{G}{2 \nu_\gamma} \qquad (0 < \nu_\gamma < d_\gamma) \\ \\ \end{align*} \]

• Applying Qi and Guo (2013) Lemma 1 to the trigamma function \(\psi^{(1)}\) shows \[ \begin{align*} \frac{d^2}{d \nu_\gamma^2} \log p(\nu_\gamma \mid \cdots) < 0 \qquad (0 < \nu_\gamma < d_\gamma) \end{align*} \]

• Similarly, \[ \begin{align*} \frac{d^2}{d \nu_\rho^2} \log p(\nu_\rho \mid \cdots) < 0 \qquad (0 < \nu_\rho < d_\rho) \end{align*} \]

• Lay out a general strategy for accelerating Markov chain Monte Carlo algorithms with parallel computation.
• Use the heterosis model as a case study.
  • Full conditionals of gene-specific parameters have general forms that can be slice-sampled in parallel.
  • Other full conditionals depend on sums and products that parallel reductions can compute quickly.
  • Carry out a proof of concept using data simulated from the model.

• Introduce an empirical Bayes counterpart method.
• Is the fully Bayesian approach really better?
  • Paschold data results
  • Heterosis gene detection power
  • False discovery rate control
  • Runtime

• The model may improve with Laplace, Student t, or horseshoe prior distributions on \(\phi_g\), \(\alpha_g\), and \(\delta_g\).
• Also borrow outside methods of analyzing heterosis data and compare as in paper 2
  • edgeR (Robinson, MD and McCarthy, DJ and Smyth, GK 2010)
  • baySeq (Hardcastle 2012)
  • Empirical Bayes for microarray (Ji, Tieming and Liu, Peng, and Nettleton, Dan 2014)
  • Empirical Bayes for RNA-seq, different model (Niemi, Jarad and Mittman, Eric and Landau, Will and Nettleton, Dan submitted 2015)